Reciprocity Theorems for Ab Initio Force Calculations 
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We present a method for calculating ab initio interatomic forces which scales quadratically with 
the size of the system and provides a physically transparent representation of the force in terms of 
the spatial variation of the electronic charge density. The method is based on a reciprocity theorem 
for evaluating an effective potential acting on a charged ion in the core of each atom. We illustrate 
the method with calculations for diatomic molecules. 



Presently, first-principles computational methods can 
be used to study the equilibrium ground-state structures 
and transient excited-state relaxation pathways even for 
relatively complex systems containing as many as sev- 
eral hundred atoms jlj. In these calculations one is 
guided through a large phase space of possible atomic 
configurations by following the gradients with respect to 
the classical nuclear coordinates of the total energy of a 
system of interacting electrons and ions. On the Born- 
Oppenheimer ground-state surface these gradients can be 
obtained by exploiting a force theorem (the Hellmann- 
Feynman theorem) which states that for any linear 
variation of a control parameter A in the quantum me- 
chanical Hamiltonian 
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where U is the total energy of the system, H is the Hamil- 
tonian and the brackets denote an expectation value in 
the electronic ground state. Using the control parame- 
ter A to denote the nuclear coordinate Ri a the theorem 
immediately provides the force acting on the I-th lattice 
site with the a-th polarization. 

In this paper we discuss a reciprocity theorem which 
provides a particularly efficient scheme for evaluating the 
forces in equation (1). One often regards (1) as measur- 
ing a response of the quantum electronic charge distribu- 
tion to an ionic displacement, and indeed the expectation 
value requires an average of a nuclear deformation poten- 
tial in the electronic ground state. Alternatively, the left 
hand side of equation (1) is a force on a classical ion re- 
sponding to the total electric field seen at the ion site. 
One can exploit this latter point of view to greatly sim- 
plify the calculation and interpretation of this force. We 
present a reciprocity theorem for inverting this problem 
which applies even in the general case where the interac- 
tion between the electronic and nuclear degrees of free- 
dom is described by effective potentials [0. Using this 
theorem the evaluation of (1) for a system containing N a 
atoms scales || as N 2 rather than as N% as in current 
electronic structure calculations Q. The method also 
greatly aids the interpretation of these forces by directly 



relating them to the spatial distribution of the ground- 
state charge density. 

Our starting point is a Hamiltonian Q describing an 
interacting system of ions and electrons, with the ionic 
degrees of freedom treated classically, H = Hk + H ee + 
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The potential Vei-km may be either a Coulomb potential 
or a more sophisticated effective potential for treating 
only the valence electronic degrees of freedom. In either 
case the Hamiltonian for the system depends parametri- 
cally on the ionic coordinates Rj only through the last 
two terms on the right hand side of equation (2) so that 
by using the force theorem, one has: 
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If we assume that the electrons and ions interact through 
the simple Coulomb potential, then one can rewrite equa- 
tion (3) from the point of view of the classical ions by first 
isolating the part of the total energy which depends on 
the nuclear coordinates: 



2 / Pion{r)V ion (r)d 6 r + 



Pci(r)V ion (r)d 3 r 
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where p c \ is the electronic charge density, Pi on {r) — 
Zre5(r — Ri) denotes the ionic charge density, and 
Mon (t) is the electrostatic potential produced by the ionic 
source distribution pi on . The second term can be rewrit- 
ten by observing that the electronic Hartree potential sat- 
isfies — j^V 2 Vh = Pd( r ): so that after integrating twice 
by parts, one has 
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C/ei-ion= / V H (r)(-— V'V ion {r))d 6 r 
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For a system of point ions, the second factor in the in- 
tegrand is Zje5(r — Ri). Thus differentiation with 
respect to the ionic coordinates explicitly shows that 
the force on each ion is a response to the electric field, 
— VVff(r) acting at the ionic site as one expects. 

The first term in expression (4) also contributes to the 
electric field acting on the 7-th site through the gradient 
— V^j_^ J Vi on (r — Rj) i.e. the gradient of the potential 
produced by the unscreened ions excluding the contribu- 
tion from the 7-th site. Combining these two contribu- 
tions for the J-th ion, one has 

V I {r) = Y,V™{r-R J ) + V H (r) (6) 

so that Fi a = —ZieVi a Vi 

The restriction on the sum in equation (6) represents 
an awkward computational constraint since in principle 
the sum is different for each ionic site in the system. How- 
ever, this can be dealt with efficiently by redistributing 
each ionic point charge uniformly on a spherical shell at 
some desired atomic sphere radius, a c in the first term 
of (6). This shell provides a potential V s which has the 
correct 1/r dependence for r > a c and is constant for 
r < a c , so that the spatially varying part of V s near the 
i-th site correctly represents the constrained sum in (6). 
With this replacement the effective potential Vj which 
governs relaxation of the classical ions can be obtained 
for all the lattice sites I in a single calculation. 

We illustrate the method with a calculation on the di- 
atomic molecule H2. In figure 1 we present a map of the 
equipotcntials of equation (6) within an atomic sphere 
centered around one of the H ions in the molecule. The 
three panels give the equipotentials calculated using the 
local density approximation (LDA) to the exchange and 
correlation energies, for three different interatomic sep- 
arations. In all cases a minimum in this potential lies 
near the center of the atomic sphere. For the compressed 
H2 the minimum lies at larger separation, and for the 
expanded molecule at a smaller separation. The force, 
calculated from the gradient of this potential, is plotted 
as a function of interatomic separation in figure 2, where 
we overlay the forces obtained by direct use of the force 
theorem (1), and as expected the two agree exactly. The 
inner panel gives the offset between the interatomic sep- 
aration and the separation which minimizes the effective 
potential computed for each separation. This shows quite 
accurately the equilibration of the molecule at the LDA 
bond length of 0.78A. 

The essence of the force theorem is that to lowest order 
in the nuclear displacements, the electronic charge den- 
sity can be regarded as rigid. In (I) one then performs an 
average of the deformation potential in this rigid density. 
In (5) we use this rigid charge density as a source term 
in the Poisson equation to extract the effective potential 
near the ionic site. The equivalence between these two re- 
ciprocal points of view applies to all types of interatomic 
interactions, including metallic, covalent, and even van 



der Waals bonding, as long as an accurate representation 
of the charge density is in hand. 

This method can be generalized to a system of elec- 
trons and ions interacting through any local pseudopo- 
tential V ps (r). A pseudopotential generally docs not re- 
tain the 1/r singularity as r — > 0, and thus the partial 
integrations leading to (5) no longer project the Hartree 
potential exactly onto a lattice site. However, by replac- 
ing Vi on by V ps in (5) one sees that the interaction en- 
ergy can always be obtained by integrating the Hartree 
potential over an effective ionic charge density which is 
calculated by taking the Laplacian of V ps . Alternatively 
we observe that the interaction energy U ps has the form: 

[/ei-ion = E / Pei(r)Vk(r ~ r 'W ~ Ri)d 3 rd 3 r' (7) 

so that by integrating over the coordinate r first, one has 

Cei-ion = /Zf y^{r')Zie5{r' - R^r' (8) 

Equation (8) explicitly shows that after the spatial av- 
erage, a point nucleus experiences an effective potential 
where pd is the source term, and for which the interac- 
tion kernel is the effective pseudopotential. If the pseu- 
dopotential is replaced by the bare Coulomb potential, 
the effective potential acting on a lattice site is just the 
electronic Hartree potential as we found earlier. 

In figure 3 we apply this method to molecular H 2 , but 
now calculated replacing the H ions by pseudo-ions. In 
the main part of the figure we overlay the forces cal- 
culated from the force theorem and calculated using the 
reciprocity relation, and again the two agree exactly. The 
inset gives the distance from the center of the sphere to 
the potential minimum as a function of interatomic sep- 
aration, and confirms that this offset crosses zero at the 
expected equilibrium separation. 

Accurate first-principles pseudopotentials frequently 
require a nonlocal representation V ps (r, r'), so that the in- 
teraction energy analogous to equation (7) has the form: 

Cei-ion = E / ^n(r)V ps (r - Rj , r' - Rj) 

n,I J 

i> n (r')d 3 rd 3 r' (9) 

and the sum is over occupied single-particle states tp n . 
Thus the interaction energy involves an integral over the 
full one-particle density matrix p(r, r') rather than simply 
the charge density alone, i.e. 

C/ cl _ ion = Trp(r,r')V ps (r,r') (10) 

Inverting this relation along the lines of equations (5) and 
(7) requires the solution of a generalized Poisson equa- 
tion with a source term which is the one-particle density 
matrix. Instead, by integrating equation (9) over a single 
intermediate coordinate s we obtain 
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This potential simplifies still further if the effective nonlo- 
cal potential is expressed as a sum of separable potentials 
as in Kleinman and Bylander's construction J5[ 

Ki-ion(r, r') = Mr ~ RiW^ c (r' - Ri) (13) 



where the <j> c (r — Ri) are a set of projection functions 
centered on the 7-th ionic site, and a c are the associ- 
ated weights. Expressing the one-electron states with 
the Fourier expansion 
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the nonlocal contribution to the interaction energy is 



U, 



-i{G-G')R! 



el-ion — ^ ^ e 
GG' nk I,c 

(k + G\(t> c )a c (4> c \k + G')c* nk (G)c nk (G') (15) 



In equation (15) the electronic ground state enters 
through the density matrix elements J2 n k c nk(G)c n k(G'), 
while the plane- wave matrix elements of the nonlocal po- 
tential appearing in the sum are normally tabulated in 
an atomic pseudopotential code. Thus the generalized 
kernel which connects the one-particle density matrix to 
the nuclear coordinate is the nonlocal potential. 

As an illustration of the method we have applied the 
prescription in equation (15) to the CI2 dimer. Here the 
Cl ions are represented by nonlocal pseudopotentials ||] 
in the Kleinman-Bylander form. In figure (4) we dis- 
play the force on each ion as a function of the inter- 
atomic separation, overlaying the forces computed from 
the Hellmann-Feynman force theorem and from the reci- 
procity relation. As expected, the two agree exactly, con- 
firming the validity of the construction. 

For a system containing N a sites, the evaluation of 
[/ei-ion (and its gradients) using the reciprocity theorem 
scales 0] as N% log 2 N a rather than as N 3 as would be 
required for a direct evaluation of equation (15). This 
occurs because the Kleinman-Bylander projection func- 
tions <j) c {r — Ri) for different atoms I, but the same atom 
type, differ only in the location of the guiding centers, Ri. 
Then one can rewrite equation (15) as a single sum 
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where the kernel f2 requires the convolution of 
4> c (G)c n k{G) with its complex conjugate, a computation 
which can be completed in 0(N a log 2 N a ) steps using the 
fast Fourier transform. Thus the evaluation of the en- 
ergy and its gradients for all sites can be accomplished in 
0(N% log 2 N a ) steps. This is a significant advance since 
it leaves the wavefunction orthogonalization operation as 
the only remaining N 3 scaling operation in plane wave 
based electronic structure calculations |l]]. 

However, in addition to its computational efficiency, a 
significant advantage of the reciprocity rule is that it re- 
lates the forces seen by the ions to the spatial distribution 
of the electronic charge density in a direct and physically 
transparent way. This is likely to be extremely useful for 
interpreting atomic relaxations near surfaces and defects 
in solids, for which the relaxations can be quite com- 
plicated, but which are nonetheless determined by the 
spatial redistribution of the valence electron density. In 
principle one might be able to exploit the composition 
of the nonlocal potential in equation (15) to clearly iso- 
late contributions to the force due to the s- p- and d-like 
components of the charge density. 

In summary, by exploiting a reciprocity relation for the 
electron-ion interaction one may reformulate the prob- 
lem of calculating atomic forces in electronic structure 
codes into a single computation which will simultane- 
ously yield the force distribution at all lattice sites in the 
system. The method scales efficiently with system size 
and directly relates the computed forces to the spatial 
distribution of the electronic charge density. 
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FIG. 1. Contour plots of the effective potential of equation 
(6) seen by an H + ion in H2 calculated at bond lengths of 
0.70A (top), 0.78A (equilibrium, middle) and 0.90A (bottom). 
The equipotentials are equally spaced at intervals of 50 meV. 
The location of the H + ion is given by the bold dot. The 
string running off each panel to the right terminates on the 
other H + ion. 
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FIG. 2. The force on an H+ ion in H2 calculated us- 
ing the full Coulomb potential plotted as a function of in- 
teratomic separation. The plot overlays the results from the 
Hellmann-Feynman force theorem (open circles) and the force 
obtained from the reciprocity relation discussed in the text 
(crosses). Here a positive force represents an attractive force. 
The vertical axis of the inner panel gives the offset between 
the interatomic separation and the distance between potential 
minima calculated for each interatomic separation. 
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FIG. 3. The force on an H + ion in H2 calculated using 
a local pseudopotential plotted as a function of interatomic 
separation. The plotting conventions are the same as those 
of figure (2). 
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FIG. 4. The force on a CI ion in CI2 calculated using a 
nonlocal pseudopotential plotted as a function of interatomic 
separation. The plotting conventions are the same as those 
of figure (2). 



4 



